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ABSTRACT 


The classical Levinson-Durbin linear prediction formulas for real valued input sequences 
are exanuned and compared to the recently proposed split- Levinson formulas. Both the 
autoregressive linear predictor model and the adaptive lattice model are used to formu- 
late the new split-Levinson algorithms. A brief introduction to the theory of symmetric 
polvnonuals is presented to form the basis of the new algorithms. Computer simulations 
are used to test and compare the computational accuracy of the new algorithms for AR 
filter coefficient estimation, parameter estimation for a moving average process, and 
spectral estimation of sinusoids in white noise. Research results indicate that the new 
algorithms reduce the number of real multiplications required for a k” order AR filter 
problem by one-half, and they are applicable to both the extended Prony method of 


spectral estimation and the estimation of moving average parameters. 
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I. INTRODUCTION 
A. OBJECTIVE 


The classical Levinson algorithm is known to provide solutions to real valued, linear 
systems involving Toeplitz structures. The computational cost for these solutions, is 
known to be O(k?), where k indicates the filter order. It has recently been proposed that 
the classical algorithm may be transformed into 2 simpler algorithms, using the theory 
of symmetric polynomials, and that either of these algorithms can be used to to solve for 
the predictor polynomial of order k at a reduced computational cost. (Ref. 1: p. 470] 

These new algorithms are termed the split-Levinson algorithms because their basis 
15 formed from the concept of symmetric polynomials. These are not new in theory, but 
the application of the process to linear prediction is a new concept. Symmetric 
polynomials are based on the Barlett Bisection Theorem [Ref. 2: pp. 1074-1076], where 
a system that possesses symmetrv about a point, such as a Toeplitz matrix, can be de- 
composed into a symmetric and an antisymmetric part. The unique point of the theory 
I$ that either part may be used to solve the problem, or a combination of both parts can 
also be used in the solution. During our research we shall only consider real data se- 
quences. 

The split-Levinson case also has a lattice structure as the classical case. However, 
as will be shown, the structure of this lattice shows little resemblance to its classical 
counterpart. A derivation of a revised split lattice structure, and its recursive algorithm 
was attempted in order to represent the split lattice in a form similar to that of the 
classical structure. Although unsuccessful, the derivation procedure is presented for 
subject matter continuity. 

Computer programs have been written to implement the new algorithms and com- 
pare them to the classical algorithms. Additionally, computer programs are included to 
apply the split-Levinson algorithm for two cases, where the computational efficiency of 
the new algorithm could be of substantial benefit. These cases include the Moving Av- 
erage (MA) problem, where the parameters of a MA model must be determined from the 
given data, and an extension of the Prony method of spectral estimation, where a least 
squares estimation of the presence of sinusoids in white noise is made from the output 


data sequence. 


This thesis compares the classical and lattice structures of the Levinson recursion 
formula given in [Ref. 3: pp. 145-167], and examines not only the formulation of the 
recursion formulas for these algorithms, but also the complexity of the computations 


and the resulting structure of each of the algorithms. 


B. THESIS ORGANIZATION 

The structure of the thesis is divided into 4 chapters, including the Introduction. In 
Chapter II we will review the classical Levinson algorithms. In the first case, the algo- 
rithm is obtained using the autocorrelation function of the input sequence, and in the 
second case, it is obtained using the forward and backward error vectors of the input 
sequence. In each case we shall establish the number of real multiplications required to 
complete a k-th order recursion of the respective algorithm. As stated, the ultimate goal 
is to establish the computational efficiency of the split-Levinson algorithm over the 
classical Levinson algorithm. Chapter III deals with the derivation of the split-Levinson 
algorithms preceded by an introductory section on symmetric polynomials. As in Chap- 
ter II, both the autocorrelation function and the lattice algorithms will be developed. 
In addition. a comparison between the computational cost of the Levinson and split- 
Levinson algorithms and an attempt to define the split lattice structure in terms similar 
to the Levinson based lattice are presented. 

In the final chapter two practical applications of the split-Levinson algorithm are 
investigated. These are: (1) the MA parameter estimation problem, and (2) the extended 
Prony method. In case (1), the Levinson recursion used to determine a predictor coeffi- 
cient vector is replaced by the split-Levinson algorithm. A comparison between the test 
coefficients and the computed coefficients is presented. In case (2), an estimation of 
sinusoids in white noise is performed. Additionally, overall conclusions of the research 


as Well as proposed topics for continued thesis research are presented. 


кә 


H. THE CLASSICAL LEVINSON ALGORITHMS 


The importance of the Levinson algorithm in linear prediction theorv is well known. 
The reason to present the algorithm in its two forms is twofold: (1) to present certain 
definitions that will be required later in the development of the split-Levinson algo- 
rithms, and (2) to detail the computational complexity of the Levinson algorithm for 
comparison purposes to the split-Levinson versions of the algorithm. In the context of 


our discussion within this thesis, we shall confine ourselves to the studv of autoregressive 


modeling problems of real sequences as in Figure 1. [Ref. 3: p. 132] 





Figure l. Autoregressive Model. 


We know from linear prediction theory the augmented normal equations given bv, 


К,а = roy (2.1) 


can optimallv be solved by the Levinson algorithm, and that this algorithm can be im- 
plemented with either the autocorrelation function or the forward and backward pre- 


diction error vectors of the input sequence [Ref. 3: pp. 152-170]. 


A. THE LEVINSON ALGORITHM 
[n order to examine the computational complexity of the Levinson recursion, it is 
necessary to formulate the recursive algorithm, and to determine the number of real 
multiplications and additions required to complete the algorithm. First, construct a 
Toeplitz matrix from the sequence s(t), of length М, defined as R, =[R,_,:0 <i, j<k], 
where the elements of the matrix are the autocorrelation lags given by [Ref 2: р. 646] 
N-1-Í 


1 


Ку =, 


s(t) s(t + 2) (0:2) 


‘= 


then, the predictor vector a can be determined as a solution to the system defined by the 


matrix equations 
СА, Зад) + Го ٣ (2.3) 


where o, 1s the prediction error norm, and is defined as 


k 
сь К+ У ay R; (2.4) 
i=] 


It is recalled from linear prediction theory, that given a positive-definite matrix R, 
of order k+1, the kth order a coefficient vector can be computed recursively from the 
nested Toeplitz submatrices, and their respective successive predictor vectors, a. The 


well-known Levinson recursion formula is this solution, and has the form 


ayy = а) + Рікі [= 012... (2:3) 


with the conditions that a,, = l, and a, ,, - 0. The parameters, p, — a,,, are called re- 


flection coefficients, also PARCOR coefficients, because they represent the partial cor- 


relation between the zero-th and the k-th element of the prediction vector with the effect 
of all the intermediate elements removed.(Ref. 4: p. 53] 
To construct the Levinson recursion we must use the prediction error norm re- 


lationship 
ср = (1 – рі)сџ (2.6) 
and the identitv 
5-1 


0-10 = = Карлу (2. 


іш0 


ty 
~ 


to define the recursion variables. Consider the following definition as it applies to the 
evinson recursion [Ref. 1: p. 472]. 


k—1 
A » Крај (2.8) 
i=0 
2 6 ر0‎ 


and solving for p, from Eq.(2.8), we have 


Ak 
Dk — 7g 





= (2.9) 
The error norm o, can be written in terms of the the normalizing term /, by rewriting 
Eq. (2.6), and making a substitution from Eq. (2.9) 


= ЦЕ б с 
= Ok-1 — یرم ےر میرم‎ (2.10) 
= o-1 кёк 


Combining Eqs. (2.5), (2.8), (2.9), and (2.10), we have the basis for the Levinson algo- 
rıthm, and it is summarized in Table 2 of Appendix A. 


B. LEVINSON LATTICE REALIZATION 

If we are given a real sequence of signal values s(0).s(1), ..., (N-1), and it is known 
that s(t) = 0, for —1 > t and ¿> N, then for the linear prediction problem of order k we 
find it necessary to find a set of real numbers a4, a4, ..., Az that will minimize the for- 


ward and backward prediction errorvectors using a linear combination of the past signal 


Un 


vectors. If we call the forward prediction error vector f(r) and the backward error vector 
bt), and define them in terms of the a, coefficients, (Ref. 2: p. 646] we have 


k 


А) = 2. s(t — i) (2.11) 


іш0 


k 
ым) = У аа 0-0 (Б 
i=0 
then it turns out that the same real numbers, a,, , will provide the solution to either of 
the forward or backward prediction problems, (ie., minimize the squared Euclidean 
norm of both /, and b,). 
Let c, be defined as the squared norm, that is 


o, = (lf, = fff, (2.13) 


From Eqs. (2.11) and (2.12) forming the first three terms of each error vector we have 


the following, 


РАО) = ayos(0) + ay,8( 1) + ... + m -Ю) (2.14) 
РМ) = а) + ay, S(0) + ajos 1) + .... + Agus 1 А) (2.15) 
(2) = аш5(2) + аз) + адоз(0) +... + а5(2 — К) (2.16) 
and, 
50) -аш(0) ға, —1) Вада (-2) +. + о ЗК) (2.17) 
ВИ) = а) + а, 1500) + ар 301) +... + аоз(1 – К) (2.18) 
602) = аш5(2) + ауу_5(1) + а, у_›5(0) +... + а5(2 — k) (2.19) 


If we examine the elements of these two vectors we can See they are related in that 
each can be derived from the other by reversing the order of the a coefficients. If we form 
the Euclidean norm of each vector, |/f,|| and |jb,!|, we see that the k-th predictor vector 
Га, 4, +---s4,.]° Minimizes the error norm, and |А| < |2,|- 

From [Ref. 3: pp. 156-157], we can use the Levinson algorithm to define the 


recursion formula for the forward and backward prediction errors given by 


А) = f(d ae Риби 11 = 1); 


bi) = Pak) + led (2.20) 


If we let the following definition apply to the lattice version of the Levinson algo- 


rıthm 


N+k-2 
= 2, ٢ (2.21) 
=1 


shencusing Eqs. (2.9), (2.10), (2.20), and (2.21) we can summarize the Levinson lattice 
algorithm in Table 3 of Appendix A. 

Even though the lattice algorithm is implemented directly from the data samples, its 
computer implementation will be more complex because of the vector manipulations 
that must occur in each iteration. The lattice structure defined by Eq. (2.20) is shown in 
Bieure 2. 

In summary, we discussed the Levinson algorithm which uses the autocorrelation 
elements ın its recursion and the related lattice structure which uses the input data di- 
rectly in its formulation. In terms of the computational complexity, both algorithms re- 
quire real multiplications of the order &?, as detailed in Tables 2 and 3 of Appendix A, 


in order to realize a K^ order predictor filter. 





Figure 2. Lattice Prediction Error Filter. 


Ш. THE SPLIT-LEVINSON ALGORITHMS 


The split-Levinson algorithms are based on the theory of symmetric and antisvm- 
metric polynomials. We know that for an k-th order real autoregressive process, the 


normal equations are 


Ко R, К, Е Ry 1 G 
P NET dn m 0 
В, К, Ro Зет Ко» а; 0 
== (3.1) 
R, ہے‎ К; 500 Ro a. J 
0 
Ra, = [o 0.0,....0J7 (3.2) 


Using the Barlett Bisection Theorem [Ref. 5: pp. 1074-1076], and because of the sym- 
metry ofthe autocorrelation matrix, we can say that the predictor coefficient vector is 


the linear combination of a symmetric and antisymmetric predictor vector given by 
а, = а + аі (3.3) 


The symmetric and antisymmetric vectors are defined as 


B 
(s) 
a = 
: гі 


3.4 
e- [5] ч 
==! 


where B represents one-half of the vector components of a,, and B' represents the re- 
versal of the vector components of B. Using Eqs. (3.3) and (3.4), we can transform the 


normal equations into 


Ra? PCM 


3.5 
Ха qe mmc DE = 


Therefore, we can see some favorable consequences of these revised normal 
equations, and their solutions. First, either the symmetric or antisymmetric form will 
give the same solution. and second, because of the symmetry of the predictor vectors, 
we need onlv solve for one-half of the predictor coefficients. 

Similar to the Levinson algorithm we now proceed to develop the split-Levinson 
algorithms from the input sequence autocorrelation function and the predictor error 


vectors. 


А. SPLIT-LEVINSON ALGORITHM 


The predictor polynomial a,(z) is defined as 


k 
-і 
ado) = У аш (3.6) 
i=0 
relative to the given Toeplitz matrix of autocorrelation lags. Denote the reverse of our 
predictor polynomial as a(z) = z-‘a,(z~'), and the predictor polynomial has been shown 


to obey the recursion (Ref. 3: 5p) 1o6- 197] 
a,(z) = a,_,(z) + ом: 1а(2) (3,7) 


and the reverse polynomial of Eq. (3.7) is 


ака) 9 z^ ay (2) - exa iz) (3.8) 

We now want to form a new polynomial from the given predictor polynomial that 
will form the basis of the split- Levinson algorithm. It is desired to show that the deter- 
mination of the coefficients of this polynomial will allow us to recover the original pre- 
dictor polynomial, and at the same time be more computationally efficient. We define 
the symmetric polynomial as P,(z), and the antisymmetric polynomial as Pz) , and we 
desire them to be of the form [Ref. 1: p. 472] 


10 


k 
P,(z) = p риб 


(а) (а) -! 
Py СЕХ Pri 2 
i=0 


Recall from Eq. (3.4) and (3.5) that the symmetric and antisymmetric predictor coeffi- 
cients are composed of two vectors that are reverses of each other, and we will define 


these vectors so that they obey the relationships 


Pri © Pk k-i 
(a) (a) 
ku Dr 


(3.10) 


Consider the mathematical interpretation of making the autocorrelation matrix, R, 
a singular matrix. If the reflection coefficient p, is made + 1, then this corresponds to 
an element of R, making the matrix singular. For this reason we shall designate the 
symmetric and antisymmetric predictor polynomials as singular predictor polynomials 
[Nek 2: p. 472] and from Eq. (3.9) thev are defined as 


=lA 


Posa ) a ez) 


PO) = a, رر‎ а 


(S: ET) 
d, (2) 


Also, these singular predictor polynomials are self-reciprocal [Ref. 2: p. 472] because of 


their symmetry and may be expressed in the following forms 


Ра) PES 


3:312 
pt? ٣٦ (3.12) 
From Eq. (3.11) we have 
=] A 
= P,(z) - 2 
PA = 2 = a 
=” — Po 12) 
If we add Eqs. (3.7) and (3.8), and make a substitution from Eq. (3.13) we have 
а(2) + 402) =2 ‘а, (2) + риа (2) + аки (2) + Рик (2) 
= (1 + ор)а, (2) + (1+ р: "4% (2) (3.14) 


= 72) 


where we have defined 7, 45 
م + 1= م2‎ (3212) 


In a similar fashion we can solve for the antisymmetric normalized singular predictor 


polynomial by subtracung Eqs. (3.7) and (3.8), and substituting from Eq. (3.13) we have 
ацг) – аг) = Ay РІ (2) (3.16) 

where 
AP -1-— p, (3.17) 


Similar to the predictor polynomial a,(z), we can define the singular predictor coef- 


ficient уесгогѕ іот ЕЧ. (3:9) А ТЕСЕ ٦٣ 


T 
Pr = [Pxo» oe ... E 





(3.18) 
ру = d PEIUS yo 1 
Since we want the split-Levinson normal equations to be of the form 
Бұқа, -(в,0,0,... 017 (3.19) 
К,а, = (0,0, ..., orl (3.20) 
then from Eq. (3.14) and (3.16) the singular predictor polynomials are 
a,(z a,(z) 
p A ) EM Е 
К 
3721 
РФ 2) „мз а,(2) ( ) 
= jm 


Since a,(z) is a polynomial formed from the predictor coefficient vector that is a solution 
to Eq. (2.3), it follows that P,(z) and Pz) are solutions to the Toeplitz system described 
bv Eqs. (3.19) and (3.20). [Ref. 1: p. 472] 


а а 
Күр, = | — | 








Ak Ak 
| (3.22) 
а а 
(a) _ k k 
Кр, Е R. -(а) а - (а) 
^k Ap 


Normalizing Eqs. (3.19) and (3.20) bv 7, and 79, the split-Levinson normal equations in 


matrix form are expressed as 


D DES НИ ТВ 





и T (5-23) 
Күр. = [re '0,0, "ого г | 
where we have defined the modified prediction error norm [Ref. 2: p. 651] 
А 
Ку m 
"Ñ (3.24) 
Т 7 
a 
k 


If we expand the matrix expressions in Eq. (3.23), the modified error norms may be ex- 


pressed as finite sums of the predictor coefficients 


k 
т = x R; Pri 


Р: (3.25) 


(а) (а) 
E » Ri Pu 


i=0 


where R, is the i-th autocorrelation element of the k-th sub matrix. 

Since the symmetric and antisymmetric polynomials are closely related, we shall 
derive only the symmetric polynomial recursion equations, and then simply present the 
results for the antisymmetric case. 

The final step in the derivation is to derive a three term recursion formula for the 
symmetric polynomial. From Eq. (3.11) and (3.14) we have the surprising result that the 
predictor polynomial a,(z) can be obtained from a linear combination of successive sin- 
olar predictor polynomials [Ref. 2: p. 472). First, form P,.,{z) from Eq. (3.11), and then 
eliminate a,(z) using Eq. (3.14) 


Ре (2) з ац(2) Я zu 


= a,{z) +27 '[4, Pz) — ay(z)] 


| Е E (3.26) 
=D Ar.) 
(1 =2 aê) = P, (2) = 2 АР) 
If we replace k by k-1 in Eq. (3.26) above we also have 
(12 Цар СЕР -: 7 B8 0) (3.27) 


We now form our recursion formula by mutiplying Eq. (3.11) by (1 — z-!), and use Eqs. 


(3.15), (3.26), and (3.27) to eliminate p, and all a, predictor polynomials. 


(г да) + (1-2 %, 10) - o1 — 2 )z ^ d, (2) 
-(1-: да (2 5р,1-: ИР,2)-а, 120 
ے‎ )1- 2 l)a,_,(z) + pi[P,(z) — ay (2) — z P,(z) + z а, (21 (3.28) 
= a, — اج‎ Pkt z py] 
- ay (z2)[(1 — z”)! — юю) 


If we now substitute for (1 --дајг-), ( - г-да, (г-) Кот Ед5. (3.27) апа (3.28), апа 
eliminate p, using Eq. (3.13), we can complete the derivation 
Ред) з z 12, P,(z) = [P (z) — z 34 PY — 29) + РИО) - РАЗА — 27) 
= 2P,(z) - Az) Pd) — > 2,1Р, (z) А РС) 
t AVPV() — P (2) =: PCE: P2) (3.29) 
Р,.1(2) = (1-27 1)Р, (2) + z 2 P,  (z2)[2, — 2J 
=(1 +z )P,(z) — ayz Р, (z) 


Taking the inverse Z transform of Eq. (3.29), we have the three term recursion formula 


for the singular predictor coefficients 
Pri = Pri F Pk k-i T %KPk=1 ki (3.30) 
where the recursion parameter a, is defined as 
Ag = ),_ [2 — Al (2.31) 


We note that r, 1s determined from P,(z) from Eq. (3.23), and therefore we conclude that 


the singular predictor polynomials can be recursively computed from Eq. (3.30). How- 


ewer, the recursiton parameter a, I$ not quite in the correct form. From Eqs. (2.10), 


(15), апд (3.31) we can alternatively compute с, from [Ref. 2: p. 473] 


Tk 


157 
k—l 





х= (3.32) 

The dual relationships for the antisymmetric split-Levinson formulas can be derived 
by following a procedure similar to the one presented above. It suffices to replace the 
quantities 4,, Te Xy Pae DY their antisymmetric duals, ie., p?, and use the following anti- 


symmetric initial conditions. [Ref. 2: p. 649] 


(a) _ 


Poo = 
(а) 
Ро + | 
(а) (3.33) 
Таса 
zs = Ro 


Recursive equations for the symmetric split-Levinson algorithm are summarized in 
Table 4 of Appendix A. Examining the entries in Table 4, we see that a full iteration 
loop of the algorithm requires approximately t real multiplications. However, because 
of the symmetry of the singular predictor coefficients, we only have to perform one-half 
of these calculations. Therefore, for a k-th order filter we need to make on the order of 
k?/2 real multiplications. The 6 function in Table 4 is used to distinguish between even 
and odd orders of the indexing variable. 

The FORTRAN program SPLIT, in Appendix B, estimates the predictor coefficients 
using the Levinson and split-Levinson algorithms. Figure 3 is a graphical comparison 
between the Known test filter coefficients of SPLIT, shown by the solid curve, and the 
filter coefficients computed by the Levinson and split-Levinson algorithms, shown by the 
dashed curve. We now undertake the derivation of the lattice form of the split-Levinson 
formulas to verify the numerical complexity of that method, and to investigate the 


symmetric and antisymmetric lattice structures compared to the Levinson lattice forms. 


В. SPLIT LATTICE ALGORITHM 

We begin the split lattice derivation by introducing the symmetric and antisymmetric 
error vectors x,, X? [Ref. 2: p. 648]. If we use previously established singularity concepts, 
and substitute + 1 for p, in Eq. (2.20) for the swmmetric and antisymmetric error vectors, 


Aand XP, respectively, we have 


COEFFICIENT NUMBER 





Figure 3. Levinson / Split Levinson Coefficient Comparison. 


0۷۶۷۶٣‏ 2 )۔ 


4 ad 
x(t) 2 f, велес 0) (3.34) 


As in the split-Levinson case, we shall proceed with the derivation of the symmetric 
split lattice, and present the antisymmetric lattice results at the end of the derivation 


with anv significant changes noted. 
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If we extend the singular polynomial concept to the singular predictor coefficients, 


we Can start with the Levinson coefficient recursion formula and substitute + | for р, 
ар = арр + Ока Iki (3.35) 
and substituting for p, 
Pri = 4р1: + ару (3.36) 


If we write Eq. (3.34) representing the time index (t) with the subscript (1) we have an 


algorithm that is more easily adapted for computers. 
Хит ар В била (3.37) 


Now, comparing Eqs. (3.36) and (3.37) we have a direct correlation between the two, 
and from the split-Levinson equation for the forward error vector, we can write the dual 


split lattice equation for x(t) 


k 


x)=) past- i) (3.38) 
i=0 


Since Eq. (3.38) is in the form of a convolution sum we can apply Z transform theory 


to see if any inferences can be made 


X,(z) = P,{z)S(z) 
А) _ P,(z) ( 
SI о 
From Eq. (3.39) we can conclude that the symmetric polynomial constitutes the Z 
transform of the transfer function formed from the error vector and input sequence z 
polynomials. Now we can use the previously derived split-Levinson algorithm, and in- 
verse transform it to obtain the lattice error vector recursion algorithm. 


Repeating the split-Levinson recursion we have 


P,,,(z) =(1 + 27) P,(z) ог Ро (2) 


Е) E (3.40) 
—-Píz-Tz Р.(2)- е Р, (2) 


Multplying Eq. (3.40) by S(z) and using Eq. (3.39) for P,(z)S(z) we have 


1 


5)2(۶ (z) = SOPLE) +z P,(z) з ау (29) 


” ” = ” = r (3.41) 
ХЕ) + X,(z) + z (2) - 0 E) 


Applying the inverse Z transform to each side of Eq. (3.41) we have the singular pre- 
dictor error vector recursion formula [Ref. 2: p. 650] 


xen ES cT ME (3.42) 


The symmetric lattice structure given by Eq. (3.42) is shown by Figure 4. [Ref. 2: 
p. 650] 

From Eq. (3.32) we know that the recursion parameter a, is defined as 7,/7,_,, and 
since a, appears in the recursion formula for the singular error vector, we need to solve 


for it. We begin with the Levinson error norm equation 


2 
вк = (1 — рај 
26, 201. l F P= pk) 


Ok 
1 + px 
2o, (1-02-20 


(3.43) 





2 7 2e, (1 — фу) 


where we have substituted v, from Eq. (3.32), and 20,_,(1 — p,) is defined as |lx,||? [Ref. 
2: p 650]; 
The initial conditions for Eq. (3.42) must be examined because most cases are trivial 


except for the case of k = 0. From Eq. (3.38) we have 
х0(4) -< Ро05(0), (3.44) 


and from [Ref. 2: p. 648] we define Py = 2. 
The antisymmetric duals are very similar to the symmetric case, and can be formed 
by replacing the symmetric variable by its antisymmetric counterpart, i.e., х?(г) for 


x(t), — p, for p,, etc. The initial conditions for the antisymmetric case are given below 


x) =0 


3.45 
х) = 8 -- зе-- 1) С 


The antisymmetric lattice structure given by the antisymmetric dual of Eq. (3.42) is 
shown in Figure 5. The split-Levinson lattice formulas given by Eqs. (3.37), (3.43), and 


(3.44) are summarized in Table 5 of Appendix A. If we examine the number of real 





Figure 4. Symmetric Lattice Structure. 
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Figure 5.  Antisvmmetric Lattice Structure. 


multiplications given in Table 4 and Table 5, then we can deduce the following conclu- 
sions. The split versions of the Levinson algorithm do produce a reduction in the com- 
plexity of the calculations when compared to the classical versions. The split-Levinson 
produces a reduction of one-half, and the Levinson produces a reduction of 3 2 (Ref. 2: 
p. 645]. The FORTRAN program SLATIS, Appendix C, implements the Levinson and 
split-Levinson lattice algorithms, and a graphical comparison between the known test 
coefficients, shown by the solid curve, the coefficients estimated by the Levinson lattice 
algorithm. shown by the dashed curve, and the coefficients estimated by the symmetric 
split-Levinson algorithm, shown by the dotted curve, is presented in Figure 6. 

The split lattice structures shown in Figure 4 on page 19, and in Figure 5 on page 
20 show that the classical lattice structure appears to be lost in the new algorithm. The 
distinct advantage of the original lattice structure is the modularity of the filter. In order 
to retrieve this appealing feature of the lattice filter we shall now proceed to derive a 
revised version of the split lattice structure, and see if it can have a form similar to the 


classical structure. 


C. SPLIT LATTICE REVISED STRUCTURE 
To begin, consider the second order classical lattice structure derived from 
Figure 2 on page 8. We can write the following equations for the first and second stage 


forward and backward prediction errors, 


fin) = sin) + pıs(n - 1) 

g (n) = pis(n) + s(n — 1) 

“ + ose - 1) 00 
gan) = piln) t+ g,(n — 1) 


Now substituting the equations for f(n) and for g,(n) into the equations for f(n) and 
gn), solving for the second stage forward and backward prediction filter errors and 


taking the Z transforms yields the transfer functions 


Е. (2) 








S(2) = |+ z^ (p, ZI 0103) sis pif (3.47) 
and. 

G,(z) = = 

SG) = р +2 Қо) +010) +2 : (3.48) 
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Figure 6. — Levinson vs. split-Levinson Coefficient Comparison. 
which obviously yield the second order predictor transfer function A,(z) and its reverse 
version А42), Where dy = 1, a, = р,. а,, = р, — PiP: = р,(1 – р,). Now forming the third 
order svmmetric polvnomial P,(z) from our second order example, we have [Ref. І: р. 


472] 


P,(z) = Aylz) + z 2 Ay(z 7!) 


3.49 
а) + (а ғау + (а +а,)г 7 + az > a 


Pevedeline a, = 1. and (a, + a) =p . then 


Вер І 


„ә 
г 

= 
— 


E 


Define the following ternis: 
ВЕ 
G=: (ре) 
therefore, 
РС)  F,(z) + CC) 


(5522) 


Equation (3.52) defines the revised symmetric split-lattice structure, and Figure 7 gives 


a graphical representation of that structure. 





Figure 7. Proposed Split Lattice Configuration 


In order to show that the preceding classical structure can be equated with the 
symmetric polynomial derived earlier, it is now necessary to form the forward and 
backward prediction error transfer function of the new split lattice structure and com- 
pare them to the Z transform of the symmetric polynomial. Thereforc, from Figure 7 


and Eq. (3.52), we can write P,(z) as 


Р) = 1+ Ку жа (K, +27) (3.53) 


Now, if we compare Eq. (3.50) to Eq.(3.53), we see that K, =p,. But, from Eq. (39) 
we know that p, = (a, + a), and substituting for a, and a, from Eq. (3.48) and (3.49) 
yields 


K, =p = pi + pX(l — pi) (3.54) 


Let us now rederive the symmetric polynomial from the revised split lattice structure. 
From Eq. (3.50) and (3.54) we have 


P(z) = 1+ (0, — py — 21032 "+2 Up рас) і 1 (3.55) 


Since we have found that the symmetric lattice can be restructured to a form similar 
to the classical lattice, the next logical step is to find a recurrence relation for the new 
lattice. Let us consider a 5% order symmetric polvnomial to determine the step-down 


recursion procedure, given by 
P,(z) = 1 + psiz2 + psgz pu раї کې‎ (3.56) 
From Eq. (3.51) we have 


Ее) + 1 + psz ` кра” 


ж E: (Ви 
6(2) = + ра? + ро 


As per the observations made in Figure 7 on page 23 and Eq. (3.51), we have the lattice 
reflection coefficient for the second stage, K, = p, Now reduce the order of F(z) to find 
the first stage reflection coefficient using the standard inverse Levinson recursion [Ref. 
4: pp. 156-157]. 


۶۰۶) 





F(z) = „2 
l-h (3.58) 
za P51 -1 
I+ K, 
therefore 
Psi : 
K, = EK (3:39) 


Now rewriting the equations for F,(z) and F,(z) using the derived reflection coefficients, 


we have 


Zu = 1 + دک‎ 
G.(z) = z! + K 

ie) nord 2 (3.60) 
Ge =e Kilkee 


and we know that 


P,(z) = F,(z) O02) 
EXEC UE a I IKE (3.61) 
تحت‎ a 


Equating the terms in Eq. (3.56) and (3.61), we have р, = K(1 + K) and pa = Ky 
Knowing the values of K, and K,, we now can form a two stage symmetric lattice similar 
at ın Figure 7 on page 23 to implement P:(z) . However, we are interested to find 
сай recursively obtain the lower orders Р.(2), P(z), etc. or the higher orders 
В 2. (2), ес. from P. (z). In an attempt to form Р,(2), уе use the standard forward 


recursion [Ref. 3: pp. 156-157] to obtain 
F(z) = F(z) + K,z `! G,(z) 


CCK CK TKK (3.62( 
KITE ET 


and 


3 (2) =2+К(1+К))2 * 
КЕ КК (3.63) 
E ТК 


Forming the symmetric polynomial P,(z) from Eqs. (3.62) and (3.63) we have 


Ра) Ва) і Сг) 
Е rm 
+ (К, + К.К + K¡KoK3)2 "+ Ky2" (3.64) 
Зе ME RON шы 
Ев сат 


Comparing Eq. (3.64) to the symmetric form of P,(z), we have 
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ре = (K, + АК, + K,K;) 


раз 1 165‏ 
نت 0 ٹب 9+4 


۹ 


1 +25 1-15 


The one-half enters into Eq. (3.65) because we Know that for even orders the polvnomial 
coefficients are svmmetric about the center element, and they must be shared in the 
matrix equations. We shall now expand Eq. (3.65) to attempt to develop a recursive 
algorithm for the sixth order coefficients from the fifth order coefficients. Expanding Eq. 
(3.65) we have 


1 + ро 
(1 + Ps2)P61 = Psi + PsiPs2 + EE ا‎ [252 + Ро - 45421 (3.66) 


where the term in brackets is an expansion of the coefficient recursion formula, Eq. 
(3.30), for p. Collecting terms we have 


1 
(1 + ps2)Pei = Psi(l + psa) + (1 + ps2) > روم + روص]روم‎ - aspan] 


"i (3.67) 
را موا موا‎ з > Pal 
Substituting for p,, from Eq. (3.30) 
Ps2 = Paa + Pal СІРІ (3.68) 


From Eq. (3.68) we can observe that the a, recursion parameter and the number 
of previous coefficients required to be known are increasing in order, and it appears that 
a simple recursive algorithm based on the above approach is not possible. Note that 
although the new lattice structure does not appear to be order recursive, we can express 
a given order symmetric or antisymmetric lattice structure in a more conventional form. 

In summary, we know that the split-Levinson algorithm is a viable replacement for 
the classical algorithm because of its computational efficiency. We have studied both 
autocorrelation and data (or lattice) realizations of the split-Levinson algorithm. An at- 
tempt to derive a recursive split lattice algorithm yielded a classical-like lattice structure, 
but it is not recursive in order. Further investigation is necessary in this direction. We 


now need to test the split-Levinson algorithms on some signal processing applications. 
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IV. APPLICATIONS OF THE SPLIT-LEVINSON ALGORITHNI 


In this chapter, we apply the spht-Levinson algorithm in (1) the MA parameter es- 
timation problem, and (2) the extended Prony method of spectral line estimation. Before 
we take up these two applications we examine the algorithm's usefulness 1f the original 


filter has coefficient symmetry, i.e., the impulse response of a linear phase FIR filter. 


А. HANKEL AND TOEPLITZ MATRICES 
іп previous derivations we have assumed the FIR filter equation to be non- 
symmetric. Let us now investigate the problem where the filter equation is symmetric, 


i.e., of the form 


y(n) = s(n) + a, s(n — 1) + “80 = 2) ЕС 


(4.1) 
۹ ۳ٰ ال‎ 
Bv definition, a symmetric polynomial is self-reciprocal, that 15 
^ =k _ y =1 
a(z) = a(z) =z a(z ) (4.2) 


Therefore, from the Levinson algorithm, predictor polynomials are known to obey the 


recurrence relation 
т 
ац2) — ay (2) + г а |(2) (4.3) 
and in our special case we have 


аа) зар) Я =a) 


a4 
-(1-рш За (2) ча 


In order to formulate a set of equations similar to the split-Levinson, 1t 15 necessary 
to derive the normal equations for our special case, and compare them to the standard 
equations, in order to develop the recursive algorithms. Since the predictor coefficients 
in a recursive algorithm produce estimates of s(n) as the algorithm steps through its re- 


cursive steps, denote this estimate as 5(п). Іп vector form, we then have 


27 


S(nln — 1) = — [s(n= 1٤ (4.5) 


To derive the normal equations we must find the minimum mean squared error (MSE) 
from the equation for the error, 


e(n) = s(n) — S(n|n — 1) (4.6) 


The minimum mean squared error is found by squaring the error term, and then differ- 
entiating the squared term with respect to the given a, vector. Combining these two ev- 
olutions we have the following equations, 
М5Е= 7 
= ETe°(n)] (4.7) 
= EL (s(n) - Qn — Y] 


To obtain the normal equations, we carry out the following steps: 


202 -2Éfs(n)s,  Te2E[S 15, 11 
ЕО 15 Да- Ел)», 0 (4.8) 


re)‏ د ار 


MIN 
E 


From the split-Levinson recursion formulas we know that the singular symmetric 


polynomial, P,(z), is defined as the following, 
Р 2) з аз е) не пра) (4.9) 


Since our predictor polynomial is symmetric, it is a reasonable question to ask if sym- 
metric polynomials also obey this recursive relationship. It is noted as an immediate 
consequence of the recursive problem, all of the preceding polynomials will also be 
symmetric. Therefore, we check the singular predictor polynomial recursion to see if it 


holds when the original polynomial is itself symmetric 
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2 u (4.10) 


=l+a2 +аг +... +2 SIS رات دہ‎ 27 


a a a ا‎ 


Now that we see that the Levinson recursive equation holds for a symmetric 
polynonual, we derive the recursive relationships for our polynomial using what is 
known from the split-Levinson equations. We have defined the symmetric polynomial 
Pz) to be a normalized combination of a non symmetric polynonual, a,(z), and its re- 


ciprocal image, a,(z) in the form of, 
Àk P, (z) = адг) + аг) (4.11) 


By direct substitution it is a trivial matter to show that this relationship also holds for 
a symmetric polynomial, a,(z). 

In order to develop the recursion for the symmetric polynomial, it is necessary to 
express the desired linear predictor in terms of the previous two predictors. To this end 


use Eq. (4.10), to form a,_,(z), and substitute from Eq. (4.11) to perform this task. 


САЛА 
а(2) ва (2)-2 а, (2) 
-1 
ayy (2) = alz) + z^ dy) 


E. (4.12) 
-а(2)%: [4,a,(z) — a,(z)] 
=a(2)(1 — z”) + z ас) 
Solving for (1 — z-!)a(z), and forming the quantity (1 — z~')a,_,(z) we have 
(1— 27 )a4G) — agi) - 2 درمز"‎ Да 
(1— z^ )aya(z) » a2) — z^ A aaa) 


These relationships will now allow us to form the three term recursion for the given 
symmetric polynomial from Eqs. (4.3), (3.11), (4.12), and (4.13) 


а(2) = а 12) + eil) 


Ж (4.14) 
а= 2٦ Закс) -аш а, (2) 


where we have defined a, 


m À, (2 — Ap) (4.15) 
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From Eq. (4.14) we can see that the coefficient recursion formula is the same form 
as Eq. (3.29), and we can deduce that the split-Levinson algorithms will work equally 
well for symmetric polynomials that describe unknown filters as it does for polynomials 


that are not symmetric. 


B. FIR MOVING AVERAGE PARAMETER ESTIMATION 
If we consider an FIR filter with an input sequence given by 


$7 — [s(n)s(n — 1)...s(n — m)], and an output y(n) given by 


M 


y(n) = ) ai s(n i) (4.16) 


=U 


then we can develop the necessary equations to estimate the moving average parameters, 
and solve for the FIR filter coefficients. The algorithm to estimate the predictor coeffi- 
cients can be defined as follows: 

Let the three predictions, y7(n), s7(m), and sf(m), represent the m-th order predictions 
of the forward estimate of y, the forward estimate of s, and the backward estimate of s 


respectively. [Ref. 6: p. 1] 


HORDIA (4.17) 
i=0 

5 (п) = c; s(n — i) (4.18) 
¿=0 

ит) = 90 - 1 )4.19( 
i=0 


where the coefficient vectors are defined as, 


bî = [1 — bo — 5. — b] 
с О NI (4.20) 


d = [0-с - 0 


Without going through all the details for the MA parmeter recursions, we can 


summarize the recursion formulas for the predictor coefficients as [Ref. 6: p. 2] 


30 


Notice that the predictor vector d" is not included in the preceding definitions since it is 
Ee reverse c” 
If we examine Eq. (4.21) we see that the recursive relationship for c" is a statement 


- 


of the Levinson recursion, ѕіпсе А”, апа А” аге the m-th order reflection coefficients. 
Therefore. we can apply the split-Levinson algorithm to solve for c^, form d" , and 
recursively determine b”. Finally, from the theory of Moving Average processes, b” = a” 
The FORTRAN program MAVI, in Appendix D, uses a 25-th order FIR test filter 
to to generate a test sequence, and the results are given in Table 6 of Appendix A. 
Figure Š is a graphical comparison between the known test coefficients, shown by 
the solid curve, and the computed filter coefficients, shown by the dashed curve. The 
curves are fitted to the linear magnitudes of the coefficients by interpolating spline 


techniques. 


C. EXTENDED PRONY METHOD 

The estimation of the existence of sinusoids in the presence of noise is a common 
occurrence in signal processing applications. In this simulation, we will show that the 
split-Levinson algorithm can be directly implemented in the process to determine the 
approximate frequencies. The noise 15 zero mean, unit variance, and white in nature. 

In this application of the split-Levinson algorithm, we attempt to approximately fit 
p exponentials to a data set of N samples. We have the constraint that V > 2p, and a 
least squares estimation procedure is used. We begin by defining the estimate of our set 
of data samples. [Ref. 7: p. 1404] 


p 
$= Dba n-0.N-1 (4.22) 


where 5, — А„ехр(}0„)/2, апа 2„=ехр(/2л/ Аг). The z, 's are roots of unit modulus with 
arbitrary frequency and occur in complex conjugate pairs as long as f,, +0 or 1/2Ar 
Therefore, in order to solve for the p sinusoids, we must solve for the roots of the fol- 


lowing polynomial. [Ref. 7: p. 1406] 
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Figure 8. MA Coefficient Comparison. 


The roots can be of unit modulus, and occur in complex conjugate conjugate pairs if 
We constrain the polynomial cocfficients to be symmetric about the center element of the 
polynomial [Ref. 7: p. 1407]. This 1s an exact ocurrence in the symmetric and anti sym- 
metric polynomials of the split-Levinson algorithm, as long as the order of the 


polynomial is even, that is 


a 
al = |: ац... 4р1 =| (4.24) 


Note that the last element of the coefficient vector is a, 2 rather than a, because of the 
symmetry of the polynomial. and that symmetric polynomials only guarantee that if a 
feet, OCCurs, then so does 1ts reciprocal z;! [Ref. 7: p. 1407]. 

The program EPRONYI, Appendix E, utilizes the split-Levinson algorithm to ap- 
proximate the p order sinusoids to the given data set. A summary of the simulation cases 


studied are presented in Table 1, and the graphical results follow. 


Table 1. UMMARY OF TEST CASES 


Ca se 


Npts 
Test hen uencies 


Ís 


Test frequencies 
Ss 
SNR 
Test frequencies 


SNR(-10, 0, 10 dB) 


fs Filter Order 
SSR 
Npts 
Test frequencies 


Л 


E SNR (-10,0,10 dB) 


SNR 
Filter Order 





1. Simulation Parameter Definitions 
All simulation cases are done in the presence of white noise, and a minimum 
possible separation frequency for the input sinusoids was determined. All plots are 
sinusoid magnitude, in a linear scale, versus digital radian frequency, 0 for0 x 0 x ٠ 
2. Simulation Results 
We begin the spectral line estimation simulations with all parameters fixed, and 
then selectively choose a parameter to vary and observe the effects of these changing 
parameters. We begin with two sinusoids of fA = 50 Hz, 4= 75 Hz, which yield 


0, = 1.396 radians, and 0, = 2.0944 radians, respectively. The number of data points 
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(`NXPTS) 1s set at 1500, and the filter order (M) is chosen го be 4 indicating the presence 
of two sinusoids. Now, we chose to vary the signal-to-noise ratio (SNR) for 10 dB. 0 dB, 
and -10 dB, shown in Figure 9. 

Figures 9 (a) and 9 (b) show that lowering the SNR from 10 dB, Figure 9 (a), 
to 0 dB. Figure 9 (b), causes the estimation to worsen, and both indicate low frequency 
estimation error. From both of these cases we can deduce the presence of two sinusoids, 
but in Figure 9 (c) it appears that only one sinusoid is present, and the low SNR has 
caused spectral estimation to fail. The conclusion for the first case is that, the better the 
SNR, the better the spectral estimate. 

For the second case we selected NPTS as the variable parameter, held the filter 
order and sinusoid frequencies constant as before, and set the SNR at 0 dB. From Fig- 
ures 10(a), (b), and (c) we clearly see that the better estimation occurs with NPTS = 
1000, Figure 10(b), because of the equal amplitude and accurate frequency estimation 
as compared to Figures 10 (a) and 10 (c). All plots show low frequency error, but also 
suggests that simulations should be conducted with NPTS set to 1000-1500 points for 
the best results. This case provides the rationale for the value of NPTS for all other 
simulations. 

In the third simulation all parameters are fixed as in the second case, and M is 
varied for 4,8, and 10. From Figures 11(a), (b), and (c) we can see that although there 
are only two sinusoids present, the estimation plots show M/2 sinusoids for all values 
of filter order. Each plot has frequency components in the vicinity of the actual fre- 
quencies, but they also give false indications of spectral lines. If we were making a deci- 
sion on the number of frequencies based on the magnitude plots of Figures 11(b) and 
(c), then signifigant errors would be introduced. In this case it is obvious that unless the 
exact number of sinusoids present 1s known, then the estimation technique will fail. 

[n the final simulation we introduce two additional sinusoids, and examine the 
effects of SNR on spectral line estimation. The constants for this case are f, = 35 Hz, 
Л - 85 Hz, f= 105 Hz, Д - 175 Hz, NPTS = 1500, M = 4, and f = 525 Hz. Digital 
frequencies are 0.419, 1.010, 1.496, and 2.094 radians per second respectively. As in case 
1, SNR takes on the values of 10 dB, 0 dB, and -10 dB. In Figure 12(a), where the SNR 
is 10 dB, we get good results with near uniform amplitude estimation, and digital fre- 
quencies that are close for all frequencies. However, from Figures 12 (b) and 12 (c), we 
can see that the spectral line for f, is missing, and an errant line appears at approximately 


2.8 radians. Both of the figures show unequal amplitude indications, but 3 of the 4 
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spectral lines are in close proxinuty to their actual values. From this case we draw the 


conclusion that the spectral line estimation performance deteriorates at low SNRs. 
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Figure 9. 
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Figure 12. Spectral Estimation (Four Sinusoids): Filter Order = 


= 8; Data Record 
Length = 1500; SNRs: (a) 10 dB, (b) 0 dB, (c)-10 dB. 
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Since the overall purpose of the simulation was to test the applicability of the 
split-Levinson algorithm to the test cases, then it has been shown that the split-Levinson 
will produce estimates for the respective cases. However, if we examine the accuracy of 
the low signal-to-noise ratio cases, we see that the new algorithm suffers a similar fate 
as the classical case in that it is difficult to accuratelv estimate the correct sinusoidal 


frequencies in the presence of noise. 


D. CONCLUSIONS 

The split-Levinson algorithm has been shown to be computationally more efficient 
than its classical counterpart. We can say that the application of the split-Levinson al- 
gorithms to practical applications in lieu of the Levinson algorithm can be advanta- 
geous, and the computational cost can be reduced significantly for large order systems. 
Additionally, the split-Levinson algorithms are applicable to problems where we want 
to model MA parameters, and perform spectral estimation using the Prony method. 

We note the following restrictive areas for the new algorithm that could make it 
unsuitable for certain signal processing applications. 

1. Non-recursive split-lattice algorithm. 


2. Computational accuracy degradation when performing spectral estimation in low 
signal-to-noise cases. 


Complexity of symmetric and antisymmetric lattice structures. 


us 


Since we have shown that the split-Levinson algorithm is a viable substitute for the 
Levinson recursion, it is reasonable to consider areas of this topic for further research. 
We know that the symmetric lattice structure can be expressed in a classical form for 
given filter orders, however, a recursive algorithm for this new structure has been elusive. 
We propose that the existing recursive algorithm for the singular predictor coefficients 
should be studied to see if a redefinition of equation parameters can extract a new re- 
cursive algorithm for the new symmetric lattice. Additionally, the algorithm's poor per- 
formance in low signal-to-noise ratio test cases of the extended Prony method is similar 
to the performance of the classical algorithm. Therefore, techniques to improve the 
classical algorithm’s performance, as detailed in [Refs. 8,9], may be investigated for ad- 
aptation to the split-Levinson algorithm. 
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APPENDIX A. TABULAR SUNIMARY OF ALGORITHMS 


The tables given in this Appendix were taken from [Ref. 2: pp.648-674], and are 


presented for the convenience of the reader. 


Table 2. THE LEVINSON ALGORITHM 


Oo d ННЯ 


For k^ 1,2,...,n 
k—1 


“Кош 2 ск 


[=0 
ао“! | рт Якір 
Oy — 0,1 — кр 
aki = Ярі Є Ребро emi 


=7 
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Table 3. THE LEVINSON LATTICE ALGORITHM 


For k= 1,2,...,n, 


N+k— 


2 
цого 2, АФ G — D 
Iz] 


Pr = ќа 
Ok = ор фу 
f) m faa) pub - 1) 


b(t) = ouf iU) t by) 
ОК) 





Table 4. THESPLIT-LEVINSON ALGORITHM 


w | |) 


3 10 
k + | = 2: — ó, 
with ó + 0 ог | 


1—1 


«= x (Ri. Ry py + дору, 
i=0 


рио = 1, 9 = трт 


Pg= d — ayl(1- Pri) 
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Table 5. THESPLIT LATIICE ALGORITHM 


x(t) = 250) Qars V=) 
x(t) = sí) + st — 1) (0 t < N) 
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Table 6. MOVING AVERAGE TEST RESULTS 
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APPENDIN B. SPLIT-LEVINSON PROGRANIS 


PROGRAM TO CALCULATE THE NTH ORDER FIR PREDICTOR FILTER USING 
THE SYMMETRIC AND ASYMMETRIC SINGULAR PREDICTOR POLYNOMIALS, 
THE SPLIT-LEVINSON RECURSION FORMULA, AND THE AUTOCORRELATION 
FUNCTION OF THE INPUT SEQUENCE: 


SIGMAN 
KEVEN 


LAMDA 


LAMDAS 


ALPHA 


ALPHAS 


KODD 


TAU 
TAUS 


RHO 
RHOS 


PS 


AS 


VARIABLE DEFINITIONS 
N-TH DEGREE NORM OF THE FILTER. 
INTEGER VARIABLE USED TO CONTROL ACCESS TO 
SUBROUTINE EVEN WHEN THE INDEX VARIABLE K IS AN 
EVEN INTEGER. 
REAL VECTOR USED WHEN DEFINING THE SINGULAR 
PREDICTOR POLYNOMIAL (PK(Z)) IN TERMS OF THE 
NORMALIZED SYMMETRIC AND ANTISYMMETRIC PARTS OF 
THE PREDICTOR POLYNOMIAL AK(Z). 
LAMDA(K) = 1 + RHO(K) 
LAMDAS(K) = 1. - RHOS(K) 

REAL VECTOR USED TO SIMPLIFY THE THREE-TERM 
RECURRENCE RELATION FOR THE SINGULAR PREDICTOR 
POLYNOMIALS OF THE SAME TYPE. 

ALPHA(K) = LAMDA(K-1)*(2. - LAMDA(K)), OR 

ALPHA(K) = TAU(K)/TAU(K-1) 

ALPHAS(K) = TAUS(K)/TAUS(K-1) 
INTEGER VARIABLE USED TO CONTROL ACCESS TO THE 
SUBROUTINE ODD WHEN THE INDEX VARIABLE K IS AN 
ODD INTEGER. 
REAL VECTOR OF "MODIFIED NORM VALUES". THE 
VALUES ARE CALCULATED FROM A SUMMATION OF 
PRODUCT TERMS OF THE AUTOCORRELATION LAGS, AND 
THE COEFFICIENTS OF THE SINGULAR PREDICTOR 
POLYNOMIALS. 
REAL VECTOR OF REFLECTION COEFFICIENTS RHO(1), 
RHO(2),... ,RHO(N). 
DESIRED ORDER OF THE PREDICTOR POLYNOMIAL. 
REAL VECTOR OF AUTOCORRELATION LAGS R(0),R(1), 
RC2Y .. 7 
ARRAY OF SINGULAR PREDICTOR POLYNOMIAL 
COEFFICIENTS FROM P(0,0),...P(N+1,N). 
ARRAY OF PREDICTOR POLYNOMIAL COEFFICIENTS. 
EACH I-TH ROW OF THE ARRAY CONTAINS THE 
PREDICTOR POLYNOMIAL COEFFICIENTS FOR THE I-TH 
ORDER PREDICTOR POLYNOMIAL. 
INTEGER VARIABLE USED IN THE SUBROUTINE ODD 
IN THE COMPUTATION OF THE TAU(K)'S. 
DUMMY VARIABLE USED DURING THE CALCULATION 
OF THE SYMMETRIC SINGULAR PREDICTOR POLYNOMIAL 
COEFFICIENTS. 


VARIABLE DECLARATIONS 


С 


С 


С 
С 


20 
С 


REAL C(0: 100),P(0: 100,0: 100) , TAU(O: 100) 

REAL A(0: 100,0: 100) , LAMDA(0: 100) , RHO( 100) , SIGMAN 
REAL PS(0:100,0: 100),AS(0: 100,0: 100) 

REAL RHOS(100),ALPHAS(100),TAUS(0: 100) , SIGMAS 
REAL AR(0:30),W(0:5000),$(0: 5000) , SIGMA(0: 100) 
REAL R(0: 100), ALPHA( 100) , LAMDAS( 0: 100) 

REAL AA(0: 100,0: 100) ,GAM( 50) , LAMK , LAM( 100) 
INTEGER M,LL,IX,T,KODD,KEVEN,L,N 


OPEN( UNIT=4 , BLANK=' ZERO’ ) 


INITIALIZE FILTER ORDER 


READ(4,100)N 
LL = N 

IX = 1913 

М 3000 


WRITE(6,300) 
WRITE(6,400)N 
WRITE(6,450)M 


DO 1 I=0,M 
CALE GAUSSCIK 1 0. V) 
ИСТ) = V 

CONTINUE 

CALL INPUT(LL,W,AR,S,M) 


INITIALIZE AUTOCORRELATION VECTOR C 


CALL ACORR(C,S,N+1,M) 
WRITE(6,2100) 


INITIAL CONDITIONS FOR THE SYMMETRIC AND ASYMMETRIC PREDICTOR 
POLYNOMIAL CALCULATIONS. 


P(0,0) = 2. 
P(1,1) = 1. 
ТАО(0) - С(0) 
LAMDA(0) = 1. 
KODD = 1 
KEVEN = 2 
A(N,0) = 1.0 
PS(0,0) = 0. 
PS(1,1) = -1. 
TAUS(0) = C(0) 

LAMDAS(0) = 1. 

AS(N,0) = 1.0 

CALL LEV(C,GAM,N,AA,LAM,SIGMA) 
WRITE(6,1700) 

WRITE(6, 1800) 

DO 20 J=1,N 
WRITE(6,1900)N,J,SIGMA(J),LAM(J),AA(N,J) 
CONTINUE 


SYMMETRIC & ASYMMETRIC SPLIT-LEVINSON RECURSION 


WRITE( 6,800) 
WRITE( 6,850) 
WRITE(6,900) 
DO 99 K=1,N 
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С CALL 


C LOOP 


P(K,0) = 1. 
РЅ(К,0) = 1. 
STATEMENTS FOR EVEN OR ODD VALUES OF INDEX K 
IF(K. EQ. KODD)THEN 
CALL AODD(C, PSK ۷٣٦ 
CALL ODD(C,P,K,N,TAU,KODD, T) 
ELSEIF(K. EQ. KEVEN) THEN 
CALL AEVEN(C,PS,K,N,TAUS, T) 
CALL EVEN(C,P,K,N,TAU,KEVEN,T) 
ENDIF 
ALPHA(K) = TAU(K)/TAU(K-1) 
ALPHAS(K) = TAUS(K)/TAUS(K-1) 
TO CALCULATE SINGULAR PREDICTOR COEFFICIENTS 
DO 40 I=1,T 
P(K+1,1) = P(K,I) + P(K,I-1) - ALPHA(K)*P(K-1,1-1) 
PS(K+1,I) = PS(K,I) * PS(K,I-1) - ALPHAS(K)*PS(K-1,I-1) 


C | DECISION PATH TO CALCULATE SYMMETRIC SINGULAR PREDICTOR COEFFICIENTS 


50 
40 


IF(I.LT.T .OR. I. EQ KF COTO 9 
L = K+1 
DO 50 J=I+1,K 
P(D.J) = P(L BT) 
PS(L,J) = -PS(L,L-J) 
CONTINUE 
LAMDA(K) = 2. - (ALPHA(K)/LAMDA(K-1)) 
LAMDAS(K) = 2. - (ALPHAS(K)/LAMDAS(K-1)) 


C REFLECTION COEFFICIENT CALCULATION 


99 


КНО(К) = LAMDA(K) - 1. 

RHOS(K) = 1. - LAMDAS(K) 
WRITE(6,1000)K,RHO(K) , RHOS(K) , GAM(K) 
CONTINUE 


C CALCULATION OF N-TH ORDER NORM (SIGMAN) AND N-TH ORDER PREDICTOR 
C COEFFICIENTS, ACN ھ۰‎ 4)٣ 


SIGMAN = LAMDA(N)*TAU(N) 

SIGMAS = LAMDAS(N)*TAUS(N) 

WRITE(6,1100) 

WRITE(6,600) 

WRITE(6,1200) 

DO 60 I=1,N 

A(N,I) = A(N,I-1) + P(N+1,1) - LAMDA(N)*P(N, I-1) 

AS(N,I) = AS(N,I-1) + PS(N+1,1) - LAMDAS(N)*PS(N,I-1) 
WRITE(6, 1300) I, TAUC I) , TAUS(I) , ALPHA(I) , ALPHAS( I) , P( I*1, I), 


%Р5(1%1,1),А(М,1),А5(М,1) 


60 

100 
200 
300 
400 
450 
600 
700 


CONTINUE 
FORMAT( I2) 

FORMAT(F12. 6) 

FORMAT('1') 

FORMAT(' FILTER ORDER = ',13) 

FORMAT('-',' NUMBER OF SAMPLE POINTS = ',15) 
FORMAT('-',103X,'FILTER COEFFICIENTS) 
FORMAT(5X,13,11X,F10. 4. 21X F104) 
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800 FORMAT('-',21X,' REFLECTION COEFFICIENTS') 
850 FORMAT(' -* ,25X,'SPLIT-LEVINSON' ) 
900 FORMAT(5X,'INDEX' ,8X,' SYMMETRIC' ,9X,' ANTISYMMETRIC' ,9X, 
+'LEVINSON' ) 
1000 FORMATS ІЗ 10Х Е12.6,12Х,Е12.6,12Х,Ғ12.6) 
1100 FORMAT('-',5X,'FILTER PARAMETERS FROM SPLIT LEVINSON RECURSION') 
1200 FORMAT(SX, INDEX 6K, TAUCK)‘ ,4X, TAUCK) AK, АБРНАСК)',4Х, 
+' ALPHA*(K)' ,4X,'P(K)',4X,'P*(K)',6X,'SYMMETRIC' ,6X, 'ASYMMETRIC' ) 
1300 FOR СИ tr l2.6,9%,Fi2. 6,ox,F12,6,4%,F 12, 6,5%,F12.6,2X, 
+F12.6,2X,F10.4,5X,F10. 4) 
1700 FORMAT('-',5X,'FILTER PARAMETERS FROM LEVINSON RECURSION' ) 
1800 FORMAT('-',8x,'INDEX' ,8X,'SIGMA(K)',5X,'LAMDACK)' ,8X,'FILTER 
+COEFFICIENTS' ) 
1900 FORMATOSK 2 ls. 7.6 12. 6,6X,F12.6,12X ,F12. 6) 
2000 FORMAT(5X, 'SPLIT-LEVINSON RECURSION CALCULATIONS’ ) 
2100 FORMAT(5X, 'INDEX' ,5X,'KNOWN COEFFICIENTS' ,5X, 'AUTOCORRELATION 
+FUNCTION C(K)') 
2200 FORMAT( 2X,13,4X,F12. 6) 
2300 FORMAT(5X,13,40X,F12. 6) 
WRITE( 6,300) 
END 
SUBROUTINE AODD(C,PS,K,N,TAUS,T) 
THIS SUBROUTINE CALCULATES THE ANTISYMMETRIC "MODIFIED 


NORMS" WHEN THE INDEX K IS AN ODD INTEGER. 
REAL C(0: 100) , PS(0: 100,0: 100) , TAUS(0: 100) 
INTEGER T 
T = (K-1)/2 
TEMP = 0. 
DO 5 I=0,T 
5 TEMP = TEMP + (C(I) - C(K-I))*PS(K,I) 


TAUS(K) = TEMP 

RETURN 

END 
SUBROUTINE ODD(C,P,K,N,TAU,KODD,T) 
THIS SUBROUTINE CALCULATES THE "MODIFIED NORMS" ( TAU(K)'S) 
WHEN THE INDEX K IS AN ODD INTEGER. 

REAL C(0: 100) ,P(0: 100,0: 100) , TAUCO: 100) 

INTEGER T 

T = (K+1)/2 

TEMP = 0. 

DO 15 I=0 T-1 

15 TEMP = TEMP + (С(Т) * C(K-ID)*P(K,I) 

TAU(K) = TEMP 

KODD = KODD + 2 

RETURN 

END 
SUBROUTINE AEVEN(C,PS,K,N,TAUS,T) 
SUBROUTINE CALCULATES THE VALUE OF THE ANTISYMMETRIC 
"MODIFIED NORMS" (TAUS(K)'S) WHEN THE INDEX K IS AN EVEN 
INTEGER. 
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са 


25 


35 


45 


15 


REAL C(O 100),PS(0: 100,0: 100) , TAUS(0: 100) 

INTEGER 7 

поко 

ТЕМР- 0. 

PS(K,T) = 0. 

DO 25 I=0,T-1 

TEMP = TEMP + (C(I) - C(K-I))*PS(K,I) 

TAUS(K) = TEMP 

RETURN 

END 
SUBROUTINE EVEN(C,P,K,N, TAU,KEVEN, T) 
SUBROUTINE CALCULATES THE VALUE OF THE "MODIFIED NORMS" 
(TAUCK)'S) WHEN THE INDEX K IS AN EVEN INTEGER. 

REAL C(0: 100) , P( 0: 100,0: 100) , TAU(O: 100) 

INTEGER T 

тив к/2 

ТЕМР- 0. 

DO 35 I=0,T-1 

TEMP = TEMP + (C(I) + C(K-IP)*P(K,I) 

TAU(K) = TEMP + C(T)*P(K,T) 

KEVEN = KEVEN + 2 

RETURN 

END 

SUBROUTINE INPUT(LL,W,AR,S,M) 

SUBROUTINE TO GENERATE THE INPUT SEQUENCE FROM A GIVEN FIR 

FILTER AND ZERO MEAN, UNIT VARIANCE WHITE NOISE. 

REAL W(0:M),AR(0: LL) ,S(0:M) 


CALCULATE INPUT SEQUENCE VALUES BETWEEN O AND FILTER ORDER. 


TEME SC 

5(0) < WCO) 

DO 45 К-1,М 

S(K) = W(K) -0.6*W(K-1) * 0.8*S(K-1) 
CONTINUE 

RETURN 

END 


SUBROUTINE ACORR(C,S,NLAG,M) 
SUBROUTINE TO CALCULATE THE AUTO CORRELATION FUNCTION OF THE 
GIVEN INPUT SEQUENCE. 

REAL C(0: 100),8(0: 5000) 
INTEGER T 

DO 5 I=0,NLAG 

TEMP = 0. 

DO 0 ٦ 

TEMP = TEMP + S(T)*S(T+1) 
CONTINUE 

CCID) = TEMP*( 1. /FLOAT ((M- 1-1) 
CONTINUE 

RETURN 

END 
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20 


30 


10 


SUBROUTINE LEV(C,GAM,N, AA, LAM, SIGMA) 

SUBROUTINE TO CALCULATE THE PREDICTOR FILTER COEFFICIENTS 
USING THE LEVINSON RECURSION. 

REAL AA(0: 100,0: 100),C(0: N*2) , GAM(N) , LAM( 100) 

REAL LAMK,SIGMA(0: 100) 

INITIAL CONDITIONS 

AA(0,0) + 1. 

SIGMA(0) + С(0) 

COMPUTE LAM(K), RHO(K); UPDATE SIGMA(K) FOR THE NEXT 
ITERATION. 


DO 10 K=1,N 

LAMK = 0. 

AA(K,0) = 1.0 

DO 20 I=0,K-1 
ГАМК = LAMK - C(K-I)*AA(K-1,1) 
LAM(K) = LAMK 

CONTINUE 


GAM(K) = LAM(K)/SIGMA(K-1) 
SIGMA(K) = SIGMA(K-1) - LAM(K)*GAM(K) 
COMPUTE THE PREDICTOR FILTER COEFFICIENTS 
01 
АА(1,1) - САМ(К) 
ELSE 
20 1-٢-٢ 
АА(К,Т) < АА(К-1,Т) + САМСК)*ААСК-1,К-Г) 
CONTINUE 
AA(K,K) = GAM(K) 
ENDIF 
CONTINUE 
RETURN 
END 
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APPENDIN C. SPLIT LATTICE ALGORITHMS 


PROGRAM TO CALCULATE THE NTH ORDER LATTICE REFLECTION 
COEFFICIENTS FROM A GIVEN SEQUENCE USING THE SYMMETRIC ERROR 
VECTOR, THE ASYMMTRIC ERROR VECTOR, OR THE FORWARD AND 
BACKWARD ERROR VECTORS. VARIABLES DEFINED IN PREVIOUS APPENDICES 
ARE NOT REDEFINED. 

vee THIS PROGRAM REQUIRES SUBROUTINE INPUT FROM APPENDIX À ww? 


SIG 
САМ 


LAM 


TAU 


AS 
AL 


AR 
XO 


LL 
X1 


AT 


TT 


tg mj c = Hr = 


VARIABLE DEFINITIONS 
N-TH DEGREE NORM OF THE FILTER. 
VECTOR OF REFLECTION COEFFICIENTS CALCULATED BY 
THE LEVINSON RECURSION. 
REAL VARIABLE USED WHEN CALCULATING THE REFLECTION 
COEFFICIENTS FROM THE LEVINSON RECURSION 
THE REFLECTION COEFFICIENT IN TERMS 
OF THE FILTER NORM IS GIVEN BY: 
RHO(K) = LAM/SIG 
REAL VECTOR OF "MODIFIED NORM VALUES". THE 
VALUES ARE CALCULATED FROM A SUMMATION OF 
PRODUCT TERMS OF THE SYMMETRIC OR ASYMMETRIC 
PREDICTION ERROR SEQUENCES. 
ARRAY OF PREDICTOR POLYNOMIAL COEFFICIENTS. 
EACH I-TH ROW OF THE ARRAY CONTAINS THE 
PREDICTOR POLYNOMIAL COEFFICIENTS FOR THE I-TH 
ORDER PREDICTOR POLYNOMIAL. 
VECTOR OF COEFFICIENTS FROM THE KNOWN TEST FILTER. 
SYMMETRIC OR ASSYMMETRIC PREDICTION ERROR VECTOR 
FOR THE (K-1) STAGE OF THE LATTICE FILTER. 
DESIRED LATTICE FILTER ORDER. 
SYMMETRIC OR ASYMMETRIC PREDICTION ERROR VECTOR 
FOR THE K-TH STAGE OF THE LATTICE FILTER. 
TEMP STORAGE FOR THE PREDICTION ERROR VECTOR WHILE 
COMPUTING THE (K+1) STAGE PREDICTION ERROR VECTOR. 
SHIFTED FORWARD PREDICTION ERROR VECTOR. 
SHIFTED BACKWARD PREDICTION EROR VECTOR. 
DESIRED ORDER OF THE PREDICTOR POLYNOMIAL. 
INTEGER VARIABLE USED IN THE PROGRAM. 
WHITE NOISE SEQUENCE VECTOR. 
INPUT SEQUENCE VECTOR 
FORWARD PREDICTION ERROR VECTOR. 
BACKWARD PREDICTION ERROR VECTOR. 


VARIABLE DECLARATIONS 


REAL AR( 30) ,W(0: 5000) ,S(0: 5000) ,RHO( 100) 
REAL A(0: 100,0: 100) ,GAM( 20) ,RHOS( 100) ,AS(0: 100,0: 100) 
REAL ALPHA,X1(0: 5000) ,X0( 0: 5000) ,AT(0: 5000) , AL(0: 100,0: 100) 
INTEGER M,LL, 1X, 0, Len 
OPEN( UNIT=4 , BLANK=' ZERO' ) 
INITIALIZE FILTER ORDER 


С) 


80 


95 


90 


READ(4,100)M 


INITIAL CONDITIONS FOR INPUT SEQUENCE GENERATOR 


LOOP 


CALL 


LL = M 
N = 5000 

Deo 
WRITE(6,200) 
WRITE(6,300)M 
WRITE(6,400)N 
WRITE(6,500) 
WRITE(6,600) 
TO GENERATE WHITE NOISE SEQUENCE AND TO READ TEST COEFFICIENTS. 
DO 1 I=0,N 

CALL ۰ ) O VW) 

МО) = У 
CONTINUE 

STATEMENT TO GENERATE INPUT SEQUENCE 
CALL INPUT(LL,W,AR,S,N) 
CALL STATEMENTS FOR SYMMETRIC, ASYMMETRIC, AND LEVINSON LATTICE 
RECURSION SUBROUTINES 

CALL SLAT(S,M,N,RHO, ALPHA ,X1, AT, XO) 
CALL ASLAT(S,M,N,RHOS, ALPHA,X1,AT, XO) 
CALL LEV(N,GAM,M,S) 
WRITE( 6,700) 
WRITE( 6,800) 
DO 80 K=1,M 
WRITE(6,900)K,RHO(K) , RHOS(K) ,GAM(K) 
CONTINUE 
DO 90 K=1,M 

ھ۰٢١‎ 

ДЕГЕ 01-1. 

АШОК 01 

IF(K .EQ. 1)THEN 

A(1,1) = RHO(K) 


AS(1,1) = RHOS(K) 
AL(1,1) = GAM(K) 
ELSE 


0 ٦۲ 
A(K,I) = A(K-1,I) + КНО(К)“А(К-Т К-Т) 


AS(K,I) = AS(K-1,1) + RHOS(K)*AS(K-1,K-1) 
AL(K,I) = AL(K-1,I) + GAM(K)*AL(K-1,K-1) 
CONTINUE 


A(K,K) = RHO(K) 
А5(К,К) = RHOS(K) 
AL(K,K) = GAM(K) 
ENDIF 
CONTINUE 
WRITE(6,1000) 
WRITECO.1100) 
WRITE( 6,1200) 
DO 96 K-1,M 
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о 


a 


WRITE(6,1300)M,K,AL(M,K),A(M,K),AS(M,K) 
96 CONTINUE 
100 ЕОКМАТ( 12) 
200 FORMAT('1') 
300 FORMAT(' FILTER ORDER = ',13) 
400 FORMAT(' ',' NUMBER OF SAMPLE POINTS = ',I5) 
500 FORMAT('-',10X,'KNOWN FILTER COEFFICIENTS!) 
600 FORMAT('-',8X,' INDEX' , 10X, ' FILTER COEFFICIENTS!) 
700 FORMAT('-',10X, REFLECTION COEFFICIENTS!) 
800 FORMAT(' -* ,5X,' INDEX ',3X,'SYMMETRIC' ,9X,'ANTISYMMETRIC' 
+,9X,' LEVINSON ") 
900 FORMAT('-',6X,I3,6X,F8. 4, 12X,F8. 4,12X, F8. 4) 
1000 FORMAT('-',15X,' FILTER COEFFICIENTS ') 
1100 FORMAT('-',20X,' LEVINSON ',12X,' SPLIT-LEVINSON ') 
1200 FORMAT(5X,' INDEX ',26X,' SYMMETRIC ',4X,' ASYMMETRIC ') 
1300 FORMAT(' ',6X,212,10X,F8. 4, 11X,F8. 4,7X,F8. 4) 
WRITE( 6,200) 
END 
SUBROUTINE SLAT(S,M,N,RHO,ALPHA,X1,AT,XO) 
SUBROUTINE TO COMPUTE THE LATTICE REFLECTION COEFFICIENTS 
USING THE SYMMETRIC PREDICTION ERROR VECTOR. 
REAL X1(0: M+N) ,XO( 0: M+N) ,RHO(M) ,S(0:N) , ALPHA 
REAL AT(0: M+N),A( 100,100) 
INITIAL CONDITIONS 
INTEGER T 
RRHO - 0. 
TAU = 0. 
INITIALIZE THE PREDICTION ERROR VECTORS FOR THE ZERO AND 1ST 
STAGES OF THE LATTICE. INITIALIZE THE ZERO STAGE MODIFIED NORM 
DO 10 T=0,N-1 
XOT) =2 sS Q) 
TAU = TAU + S(T)**2 
10 CONTINUE 
DO 20 T=0,N 
IF(T.EQ.N)S(T) = 0. 
IF(T. EQ. O) THEN 
X1(T) = S(T) 
ELSE 
Х1(Т) - S(T) * S(T-1) 
ENDIF 
20 CONTINUE 
LOOP TO COMPUTE THE REFLECTION COEFFICIENTS 
DO 101 K=1,M 
TTAU = TAU 
STORE TAU(K-1), AND COMPUTE TAU(K). 
ТАП = 0. 
DO 30 T=0,N+K-1 
TAU = TAU Е KI(T)%%2 
30 CONTINUE 
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TAU = TAU/2. 


C COMPUTE ALPHA(K), RHO(K); STORE TAU(K) AND RHO(K) FOR 
C NEXT ITERATION. 

ALPHA = TAU/TTAU 

TTAU = TAU 


RHO(K) = 1. - (ALPHA/(1. + RRHO)) 
RRHO = RHO(K) 
C LOOP 180 0011015 THE (K+1) PREDIGTION ERROR VECTOR, 
DO 40 T=0,N+K 
7 0 ۳ 
АТИ) = ХИТ) 
ظ2‎ ۶۹۳ 
В = ТЕ 


ELSE 
AED SKID сре (ТЕТ = ALPHA*XO(T=1) 
ENDIF 
40 CONTINUE 
C LOOPS TO UPDATE PREDICTION ERROR VECTORS FOR NEXT ITERATION. 
C (SHIFT X1 INTO XO AND AT INTO X1) 


DO 50 T=0,N+K-1 
XO(T) = X1(T) 
50 CONTINUE 
DO 60 T=0,N+K 
X1(T) = AT(T) 
60 CONTINUE 
101 CONTINUE 
RETURN 
END 
SUBROUTINE ASLAT(S,M,N,RHOS, ALPHA ,X1,AT, XO) 
С SUBROUTINE TO COMPUTE THE LATTICE REFLECTION COEFFICIENTS 
З USING THE ASYMMETRIC PREDICTION ERROR VECTOR. 
REAL X1(0: M+N),XO(0: M+N),RHOS(M),S(O0:N),ALPHA 
REAL AT(O: M+N) 
INTEGER T 
С INITIAL CONDITIONS 
RRHO = 0. 
TAU = 0. 
INITIALIZE THE PREDICTION ERROR VECTORS FOR THE ZERO AND 1ST 
STAGES OF THE LATTICE. INITIALIZE THE ZERO STAGE MODIFIED NORM 
DO 10 T=0,N-1 
XO(T) = 0. 
TAU = TAU + S(T)**2 
10 CONTINUE 
DO 20 T=0,N 
IF(T.EQ.N)X1(T) = -S(T-1) 
IF( T. EQ. O) THEN 


а о 


X1(T) = S(T) 

ELSE 

КЕ) = SCT) - SCT=1) 
ENDIF 
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20 CONTINUE 
C LOOP TO COMPUTE THE REFLECTION COEFFICIENTS 
DO 101 К-1,М 
C STORE TAU(K-1), AND COMPUTE TAU(K). 
TTAU = TAU 
TAU = 0. 
DO 30 T=0,N+K-1 
TAU = TAU + X1(T)**2 
30 CONTINUE 
TAU = TAU/2. 
С COMPUTE ALPHA(K), RHO(K); STORE TAU(K) AND RHO(K) FOR 
C NEXT ITERATION. 
ALPHA = TAU/TTAU 
TTAU = TAU 
RHOS(K) = (ALPHA/(1. - RRHO)) - 1. 
RRHO = RHOS(K) 
C LOOP TO COMPUTE ТНЕ (К%1) PREDICTION ERROR VECTOR. 
DO 40 T=0,N+K 
IF(T .EQ. 0)THEN 
AT(T) = X1(T) 
ELSEIF(T. EQ. №+К)ТНЕМ 
AT(T) = X1(T-1) 


ELSE 
AT(T) = Х1(Т) + X1(T-1) - ALPHA*XO(T-1) 
ENDIF 
40 CONTINUE 
С LOOPS TO UPDATE PREDICTION ERROR VECTORS FOR NEXT ITERATION. 
C (SHIFT X1 INTO XO AND AT INTO X1) 


DO 50 T=0,N+K-1 
XO(T) = X1(T) 
50 CONTINUE 
DO 60 T=0,N+K 
X1(T) = AT(T) 
60 CONTINUE 
101 CONTINUE 
RETURN 
END 
SUBROUTINE LEV(N,GAM,M,S) 
C SUBROUTINE TO COMPUTE THE REFLECTION COEFFICIENTS FOR AN 
N-TH ORDER LATTICE FILTER FROM THE FORWARD AND BACKWARD 
C PREDICTION ERROR VECTORS. 
REAL F(0:5100),BC0:5100),; FTCOS5 100) 9BTCD: 5189004 GAMC2O0) 
REAL LAM,SIG,S(0: N) 
INTEGER T 
C INITIAL CONDITIONS; INITIALIZE FORWARD AND BACKWARD PREDICTION 
C ERROR VECTORS. 
SIG - O. 
DO 10 I=0,N-1 
F(I) =S(I) 
BCI) = S(T) 


о 
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ao 


10 


30 


40 


50 


60 
20 


)۳ئ 
SIG SIG + 5СІ “2‏ 
CONTINUE‏ 


LOOP TO COMPUTE THE REFLECTION COEFFICIENTS 


DO 20 K=1,M 


FORM SHIFTED ERROR VECTORS 


DO 30 T=1,N+K-1 
BT(T) = B(T-1) 

CONTINUE 

ВТ(0) = 0. 

FT(N+K-1) = 0. 
LAY = 0. 


COMPUTE LAM(K), GAM(K); UPDATE K-TH ERROR NORM AND 
STORE FOR NEXT ITERATION. 


DO 40 T=1,N+K-2 
LAM = LAM - FT(T)*BT(T) 

CONTINUE 

GAM(K) = LAM/SIG 

IF(K .EQ. M)GOTO 20 

SIG = SIG - LAM*GAM(K) 


COMPUTE (K+1) FORWARD AND BACKWARD PREDICTION ERRORS AND SHIFT 
INTO TEMPORARY VECTORS FOR NEXT ITERATION. 


DO 50 T=0,N+K-1 

F(T) = FI(T) + GAM(K)*BT(T) 
B(T) = САМ(К)*ЕТСТ) + ВТ(Т) 
CONTINUE 
DO 60 T=0,N+K-1 

FI(T) = F(T) 

BT(T) = B(T) 

CONTINUE 

CONTINUE 

RETURN 

END 
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APPENDIX D. MA PREDICTOR COEFFICIENT PROGRAM 


THIS PROGRAM IS TO COMPUTE THE FIR FILTER COEFFICIENTS 
USING THE SPLIT-LEVINSON ALGORITHM, AND AN AUTOREGRESSIVE 
MOVING AVERAGE PROCESS. VARIABLES DEFINED IN PREVIOUS 
APPENDICES ARE NOT REDEFINED. 
‚rt THIS PROGRAM REQUIRES THE SUBROUTINES ODD, AND EVEN 
FROM APPENDIX A er 


ENORM 


NSTRT 
NEND 
DELX 
DELY 
NPTS 
EXB 
RXY 
AA 


Y 


VARIABLE DEFINITIONS 
PREDICTOR COEFFICIENT ERROR NORM. 


ENORM - 50КТ((А(0)-АА(0))%%2 +. ..+ CACN) -AA(N) )**2) 


NUMBER OF POINTS OF INPUT SEQUENCE TO START. 
NUMBER OF POINTS OF INPUT SEQUENCE AT END OF PROGRAM. 
ERROR VECTOR. 

ERROR VECTOR. 

NUMBER OF INPUT DATA POINTS (0,1,...,NPTS). 
BACKWARD PREDICTION ERROR POWER OF X. 

VECTOR OF X AND Y CROSSCORRLATION ELEMENTS 
VECTOR OF CALCULATED PREDICTOR COEFFICIENTS. 
INDEX FOR X-AXIS OF PLOTTING ROUTINE. 

VECTOR OF PREDICTOR COEFFICIENT NORMS. 

FORWARD PREDICTION ERROR POWER OF X. 

FORWARD PREDICTION ERROR POWER OF Y. 

Y REFLECTION COEFFICIENT. 

X REFLECTION COEFFICIENT. 

FILTER ORDER VARIABLE USED IN SUBROUTINE CORR. 
INTEGER SEED NUMBER USED BY IMSL SUBROUTINE GAUSS. 
X AUTOCORRELATION VECTOR. 

Y AUTOCORRELATION VECTOR. 

MA COEFFICIENT VECTOR. 

MA COEFFICIENT VECTOR. 

MA COEFFICIENT VECTOR. 

INTEGER VARIABLE USED IN THE SUBROUTINE ODD 

IN THE COMPUTATION OF THE TAU(K)'S. 

DUMMY VARIABLE USED DURING THE CALCULATION 

OF THE SYMMETRIC SINGULAR PREDICTOR POLYNOMIAL 
COEFFICIENTS. 

INPUT WHITE NOISE VECTOR. 

INDEXING VARIABLE USED IN FIR FILTER COEFFICIENT 
RECURSION (M=1,2,...,N). 

OUTPUT SEQUENCE FROM FIR TEST FILTER. 


VARIABLE DECLARATIONS 
REAL P(0: 100,0: 100) , TAU(0: 100) ,C(0: 50) ,B(0: 50), EN(200) 
REAL A(0: 100,0: 100) , LAMDA(0: 100) ,X(0: 10000) ,D(0: 50) 
REAL AR(O: 30),Y(0: 10000) ,EY(0: 50) ,EX(0: 50) ,KX(50) 
REAL DELX(0: 50),DELY(0: 50) ,EXB(0: 50) ,KY( 50) ,XX(200) 
REAL RX(0: 100) , ALPHA( 100) , RY(0: 2) ,RXY(0: 100) , AA(0: 100) 
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INTEGER M,LL,IX,T,KODD,KEVEN,L,N,NPTS 
C DESIRED FILTER ORDER AND THE TEST FILTER COEFFICIENTS 
C ARE READ FROM A DATA FILE (FILE FTO4F001). 
OPENCUNIT=4,BLANK='ZERO” ) 
C INITIALIZE FILTER ORDER AND TEST FILTER COEFFICIENTS 
READ(4,100)N 
= 
DO 20 K=0,LL 
READ(4,700)ARCK) 
20 CONTINUE 
NPTS = 4000 
WRITE(6,300) 
WRITE(6,400)N 
IX = 1913 
WRITE(6,600)NPTS 
C GENERATE (NPTS+1) WHITE NOISE SEQUENCE 
DO 10 I=0,NPTS 
CAUCA LSSI LOS, V) 
NY 
10 CONTINUE 
C CREATE INPUT SEQUENCE FROM GIVEN FIR TEST FILTER 
CALL INPUT(LL,X,AR,Y,NPTS) 
C INITIALIZE AUTOCORRELATION VECTOR RX,RY, AND CROSSCORRELATION 
C VECTOR RXY 
CALL CORR(N+1,NPTS,X,Y,RX,RY,RXY) 
WRITE(6,800) 
DO 30 I=0,N 
WRITE(6,900)I,AR(I),RX(I) 
30 CONTINUE 
C INITIAL CONDITIONS FOR MOVING AVERAGE MODEL PARAMETERS 
BOO = -RXY(0)/RX(0) 


EY(0) = RY(0) - BOO*RXY(0) 
EX(0) = RX(0) 
DO 40 I=0,1 
C(I) = I 
D(I) = I 
40 CONTINUE 
رح یج‎ 
B(1) = BOO 


DELY(0) = RXY(1) - BOO*RX(1) 
DELX(0) = RX(1) 
C LOOP TO GENERATE PREDICTOR COEFFICIENT VECTOR 
DO 120 M=1,N 
C INITIAL CONDITIONS FOR THE SYMMETRIC PREDICTOR POLYNOMIAL 
C | CALCULATIONS. 
P950) 
PCY, 1) 1. 
TAU(0) RX(0) 
LAMDA(0) = 1. 
KODD = 1 


2 
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КЕУЕМ - 2 
A(M,0) = 1. 
C SYMMETRIC SPLIT-LEVINSON RECURSION 
DO 70 K=1,M 
P(K,0) = 1. 
C CALL STATEMENTS FOR EVEN OR ODD VALUES OF INDEX K 
IF(K. EQ. KODD)THEN 
CALL ODD(RX,P,K,N,TAU,KODD,T) 
ELSEIF(K. EQ. KEVEN) THEN 
CALL EVEN(RX,P,K,N, TAU, KEVEN, T) 
ENDIF 
ALPHA(K) = TAU(K)/TAU(K-1) 
C LOOP TO CALCULATE SINGULAR PREDICTOR COEFFICIENTS 
DO 60 I=1,T 
P(K+1,I) = P(K,I) + P(K,I-1) - ALPHA(K)*P(K-1,1-1) 
C DECISION PATH TO CALCULATE SYMMETRIC SINGULAR PREDICTOR COEFFICIENTS 
ТЕСТ. БЕТ DR. 1. EO KCOTO CO 
L = K+1 
DO 50 J=I+1,K 
PCL: = EL L-I) 
50 CONTINUE 
60 CONTINUE 
LAMDA(K) = 2. - (ALPHA(K)/LAMDA(K-1)) 
70 CONTINUE 
C CALCULATION OF N-TH ORDER PREDICTOR COEFFICIENTS, A(K,1),...A(K,K) 
С COMPUTE PREDICTOR COEFFICIENTS FOR K-TH ITERATION 
DO 80 I=1,M 
A(M,I) = A(M,I-1) + P(M+1,I) - LAMDA(M)#P(M,I-1) 
С(І) - А(М,1) 


80 CONTINUE 
DO 90 J=1,M 
DCI) = -G(9) 


700 . LT. 71 
DO 95 I=J+1,M 
D(I) = D(I-1) 
95 CONTINUE 
ENDIF 
90 CONTINUE 
D(0) = O. 
D(M+1) = 1. 
@ UPDATE PREDICTION ERRORS 
XBTMP = 0. 
XTMP = 0. 
DO 25 I=1,M 
XBTMP = XBTMP + C(I)*RXY(M+1-1) 
ХТМР з XTMP * C(I)*RX(M-*1-I) 
25 CONTINUE 
DELX(M) = RX(M+1) - XTMP 
EXB(M) = RXY(M+1) - XBTMP 
С UPDATE REFLECTION COEFFICIENTS 
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KX(M) =DEEACM=1) /EXCH=1) 


EX(M) = EX(M-1) + KX(M)*DELX(M-1) 
KY(M) = -DELY(M-1)/EX(M) 
Ba) EYNI) + KY(M)ZBXB(M) 


C UPDATE B VECTOR 


45 
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150 


120 
100 
200 
300 
400 
500 
600 
700 
800 


900 
1000 
1100 
1200 


a 


B(M+1) = 0. 
DO 45 I=0,M+1 
В) = ВСІ) + KYCM)*DCT) 
CONTINUE 
YTMP = 0. 
DO 55 I=1,M+1 
УТ = YDUPBESSBODUSRX(MET2-I) 
CONTINUE 
DEGY(M) = RAYCM+1) = YIMP 
0100 НЕМ 
WRITE(6,1000)N 
WRITE(6,1100) 
ВО 130 К-0,М 
НЕТТЕ( 6. 1200)К, -В(К+1) 
AA(K) = -B(K+1) 


CONTINUE 
ENDIF 

CONTINUE 

FORMAT( 12) 

FORMAT(' ') 

FORMAT( '1') 

FORMAT(' FILTER ORDER = ',13) 

FORMAT( ' -' ) 

FORMAT('-',' NUMBER OF SAMPLE POINTS = ',15) 
FORMAT(F8. 4) 

FORMAT('-','5X,' INDEX' ,5X, KNOWN COEFFICIENTS! ,5X, 


+'AUTOCORRELATION FUNCTION R(K)') 


FORMATO (se 5Х,13,11Х,Е8 4,215 Е8.4) 
FORMAT('-',10X,'PREDICTOR COEFFICIENTS FOR FILTER ORDER = ',1I3) 
FORMAT('-',5X,'INDEX',12X,'FIR PREDICTOR COEFFICIENTS!) 
FORMATC 3 45%,13,23X,F8.4) 

WRITE( 6,300) 
END 

SUBROUTINE INPUT(LL,X,AR,Y,NPTS) 

SUBROUTINE TO GENERATE THE INPUT SEQUENCE FROM A GIVEN FIR 

FILTER AND ZERO MEAN, UNIT VARIANCE WHITE NOISE. 

REAL X(0: NPTS) ,AR(O: LL) , Y(0: NPTS) 


C CALCULATE INPUT SEQUENCE VALUES BETWEEN O AND FILTER ORDER. 


DO 45 К-0,МРТ5 
ТЕМР = 0. 
IF(KCLE. LL)THEN 
LU=K 
ELSE 
LU=LL 
ENDIF 
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ех 
DO 55 I=0,LU 
TEMP = TEMP + AR(I)*X(J) 
J = J-1 
55 CONTINUE 
Y(K) = TEMP 
45 CONTINUE 
RETURN 
END 
SUBROUTINE CORR(NLAG,NPTS,X,Y,RX,RY,RXY) 
SUBROUTINE TO CALCULATE THE AUTOCORRELATION FUNCTION X AND Y AND 
THE CROSSCORRELATION FUNCTION BETWEEN X AND Y 
REAL RX(0: NLAG) ,Y(0: NPTS) ,X(0: NPTS) , RXY(0: NLAG) , RY(0: 2) 
INTEGER T 
COMPUTE THE AUTOCORRELATION OF X AND THE CROSSCORRELATION OF 
X AND Y FOR LAGS Оо. NLE 
DO 5 I=0,NLAG 
XTEMP = 0. 
XYTEMP = 0. 
COMPUTE THE AUTOCORRELATION OF X AND THE CROSSCORRELATION OF 
X & Y FOR LAG I 
DO 15 T=0,NPTS-1-1 
ХТЕМР = ХТЕМР + Х(Т)*Х(Т+І) 
ХҮТЕМР = ХҮТЕМР + Х(Т)*Ү(Т+І) 
15 CONTINUE 
RX(I) = XTEMP*(1. /FLOAT(NPTS-1-1)) 
ВХУ(І) є XYTEMP*(1. /FLOAT(NPTS-1-1)) 
COMPUTE THE ZERO LAG AUTOCORRELATION FUNCTION OF Y 
IF(I .EQ. O)THEN 
RY(0) = 0. 
DO 16 J=0,NPTS-1 
RY(0) = ВУ(0) + Y(J)**2 
16 CONTINUE 
ВУ(0) = RY(0)*(1. /FLOAT(NPTS-1)) 
ENDIF 
5 CONTINUE 
RETURN 
END 
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ته تک 0 ال6۱ لن د رت СО (е cy‏ 


APPENDIX E. EXTENDED PRONY PROGRAM 


PROGRAM TO CALCULATE THE NTH ORDER LATTICE REFLECTION 
COEFFICIENTS FROM A GIVEN SEQUENCE USING THE SYMMETRIC ERROR 
VECTOR, THE ASYMMTRIC ERROR VECTOR, OR THE FORWARD AND 
BACKWARD ERROR VECTORS. 

VARIABLE DEFINITIONS 


РТ - TEMPORARY ARRAY USED TO AVERAGE PREDICTOR 
COEFFICIENTS. 

PP - ESTIMATED NUMBER OF COMPLEX SINUSOIDS PRESENT. 

Al = AMPLITUDE OF #1 SINUSOID, (1-4) SINUSOIDS 
PRESENT: 

FS - SAMPLING FREQUENCY. 

ЕТ = FREQUENCY OF #1 SINUSOID IN TEST SEQUENCE. 


ТНЕТА1 - DIGITAL FREQUENCY OF #1 TEST ANALOG FREQUENCY. 


VARIABLE DECLARATIONS 


REAL W(0: 5000) ,S(0: 5000) , ALPHA( 100) , ROOTR( 100) , XCOF( 0: 100) 

REAL P(0: 100,0: 100) , ALPHAS( 100) , COF(0: 100) , ROOTI( 100) 

REAL X1(0: 5000) ,X0(0: 5000) ,AT(0: 5000) ,PS(0: 100,0: 100) ,РТ(0: 100) 
INTEGER T,PP 


OPEN(UNIT=4,BLANK='ZERO' ) 


INITIALIZE FILTER ORDER 


READ(4,100)PP 


M = 2*PP 
INITIAL CONDITIONS FOR INPUT SEQUENCE GENERATOR 

IX = 1913 

A = SQRT(2. ) 
B = SQRT(2. ) 
C = SQRT(10. ) 
D = SORT(2.0) 
E = SQRT(2. 0) 
Е1 = 5.5Е+01 
Е2 = 7.5Е+02 
ЕЗ = 1.25Е+02 
Е& = 1.75Е+02 


LOOP 


FS = 2.25E+02 

TWOPI = 6.2831853 

THETA1 = (TWOPI*F1)/FS 

THETA2 = (TWOPI*F2)/FS 

TO GENERATE WHITE NOISE SEQUENCE AND TO READ TEST COEFFICIENT 
DO 1 I=0,N 

CALL GAUSS(IX,1. ,0. ,V) 

ИСТ) = C*V 


Al = A*COS( TWOPI*(F1/FS)*FLOAT(I)) 
A2 - B*COS(TWOPI*(F2/FS)*FLOAT(I)) 
АЗ 2 D*COS(TWOPI*(F3/FS)*FLOAT(I)) 
AS 2 E*COS(TWOPI*(F4/FS)*FLOAT(I)) 


S(I) = Al + A2 + W( 1) 
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CONTINUE 


C CALL STATEMENTS FOR SYMMETRIC, ASYMMETRIC, AND LEVINSON LATTICE 
C RECURSION SUBROUTINES 


C 


e 


C) C) 


11 

100 
200 
300 
400 
500 
600 
700 


CALL SLAT(S,M,N,P, ALPHA, X1, AT, XO) 
CALL ASLAT(S,M,N,PS, ALPHAS ,X1,AT, XO) 


WRITE( 6,200) 


WRITE(6,300)PP 


WRITE(6,400)M 
WRITE(6,500)N 


DISPLAY COEFFICIENTS OF POLYNOMIAL 
WRITE(6,600) 

DO 11 K=0,M 

IF(K ۹9 TOBE KRI 0 
ИКТТЕ( 6, 700)М,К, РМ К) 


CONTINUE 
FORMAT( I4) 

FORMAT( '1') 

FORMAT(' NUMBER OF COMPLEX EXPONENTIALS IN MODEL = ',I3) 
FORMAT(' ',' SYMMETRIC FILTER ORDER = ',13) 


FORMAT(' ',' NUMBER OF SAMPLE POINTS = ',I5) 
FORMAT('-',8X,' INDEX' , 13X, ' COEFFICIENTS ) 
FORMAT( 5X,212,16X,F8. 4) 


WRITE(6,200) 


END 


SUBROUTINE SLAT(S,M,N,P,ALPHA,X1,AT, XO) 


SUBROUTINE TO COMPUTE THE SYMMETRIC PREDICTOR COEFFICIENTS 
USING THE SYMMETRIC PREDICTION ERROR VECTOR. 
REAL X1(0: MN) ,XO( 0: M*N) , ALPHA(M) , S(0: №) 
REAL AT(0:M+N),P(0: 100,0: 100) 

INTEGER T 

INITIAL CONDITIONS 

KODD = 1 

KEVEN = 2 

TAU = 0. 

.1 ت۶ 

P(0,0) = 2. 


INITIALIZE THE PREDICTION ERROR VECTORS FOR THE ZERO AND 1ST 
STAGES OF THE LATTICE. INITIALIZE THE ZERO STAGE MODIFIED NORM 


10 


DO 10 T=0,N-1 

XO(T) = 2.*S(T) 

TAU = TAU + S(T)**2 

CONTINUE 

DO 20 T=0,N 
IF(T. EQ. N)S(T) = 0. 
IF(T. EQ. O)THEN 


АИТ) = SCT) 
ELSE 

IT) = SCI) = 500-0) 
ENDIF 
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20 


30 


50 
40 


60 


С) 


710, 


80 


CONTINUE 
LOOP TO COMPUTE THE SYMMETRIC PREDICTOR COEFFICIENTS 
DO 101 K=1,M 
P(K,0) = 1.0 
0801 = ТАП 
IF(K .EQ. KODD)THEN 
LL - (K*1)/2 
ELSE 
LL = K/2 
ENDIF 
STORE TAU(K-1), AND COMPUTE TAU(K). 
TAU = 0. 
DO 30 T=0,N+K-1 
TAU = TAU + X1(T)**2 
CONTINUE 
TAU = TAU/2. 
COMPUTE ALPHA(K); STORE TAU(K) FOR NEXT ITERATION. 
ALPHA(K) = TAU/TTAU 
TTAU = TAU 


LOOP TO CALCULATE SINGULAR PREDICTOR COEFFICIENTS 


DO 40 I=1,LL 
ВВ, RIE PK IS) АБРНА(К)ЗРСК-Л Те 


DECISION PATH TO CALCULATE SYMMETRIC SINGULAR PREDICTOR COEFFICIENTS 


TEM LE ОВ. 0 
L = К+1 
DO 50 J=I+1,K 
ВОРЮУ) = P(L,L-J) 
CONTINUE 
CONTINUE 
LOOP TO COMPUTE THE (K+1) PREDICTION ERROR VECTOR. 
DO 60 T=0,N+K 
IF(T .EQ. 0)THEN 
AT(T) = X1(T) 
ELSEIF(T. EQ. N+K)THEN 
АТ(Т) = Х1(т-1) 
ELSE 
AT(T) = X1(T) + X1(T-1) - ALPHA(K)*XO(T-1) 
ENDIF 
CONTINUE 
LOOPS TO UPDATE PREDICTION ERROR VECTORS FOR NEXT ITERATION. 
(SHIFT X1 INTO XO AND AT INTO X1) 
DO 70 T=0,N+K-1 
XO(T) = X1(T) 
CONTINUE 
DO 80 T=0,N+K 
X1(T) = AT(T) 
CONTINUE 
IF(K .EQ. KODD)THEN 
KODD = KODD + 2 
ELSE 
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aa 


KEVEN = KEVEN + 2 
ENDIF 
101 CONTINUE 
RETURN 
END 
SUBROUTINE ASLAT(S,M,N,PS,ALPHAS,X1,AT,XO) 
SUBROUTINE TO COMPUTE THE LATTICE REFLECTION COEFFICIENTS 
USING THE ASYMMETRIC PREDICTION ERROR VECTOR. 
REAL X1(0: M*N) ,XO( 0: M-N) , PS(0: 100,0: 100) ,S(0: N) , ALPHAS(M) 
REAL AT(0:M*N) 
INTEGER T 
INITIAL CONDITIONS 
Р5(0,0) = 2. 
PS(1,1) = 1. 
KODD = 1 
KEVEN = 2 
TAU = 0: 
INITIALIZE THE PREDICTION ERROR VECTORS FOR THE ZERO AND 1ST 
STAGES OF THE LATTICE. INITIALIZE THE ZERO STAGE MODIFIED NORM 
DO 10 T=0,N-1 
XO(T) = 0. 
TAU = TAU + S(T)**2 
10 CONTINUE 
DO 20 T=0,N 
IF(T. EQ.N)X1(T) = -S(T-1) 
ТЕСТ. EQ. 0) THEN 
X1(T) = S(T) 
ELSE 
X1CT) та (павета) 
ENDIF 
20 CONTINUE 
LOOP TO COMPUTE THE REFLECTION COEFFICIENTS 
DO 101 К-1,М 
PS(K,0) = 1. 
TTAU = TAU 
IF(K .EQ. KODD)THEN 
LL = (K-1)/2 
ELSE 
LL = K/2 
ENDIF 
STORE TAU(K-1), AND COMPUTE TAU(K). 
TAU = 0. 
DO 30 T=0,N+K-1 
TAU = TAU + X1(T)**2 
30 CONTINUE 
TAU = TAU/2. 
COMPUTE ALPHA(K) ; STORE TAU(K) FOR NEXT ITERATION. 
ALPHAS(K) = TAU/TTAU 
TTAU = TAU 
LOOP TO CALCULATE SINGULAR PREDICTOR COEFFICIENTS 


DO 40 I=1,LL 
PS(K+1,1) = PS(K,I) + PS(K,I-1) - ALPHAS(K)*PS(K-1,I-1) 
C DECISION PATH TO CALCULATE SYMMETRIC SINGULAR PREDICTOR COEFFICIENTS 
IF(I. LT. LL .OR. I.EQ.K)GOTO 40 
L = 2 7 
DO 50 J=I+1,K 
PS(L,J) = -PS(L,L-J) 


50 CONTINUE 
40 CONTINUE 
C 1006 18 081008 THE (Ktl) “PREDICTION ERROR VECTOR. 


DO 60 T=0,N+K 
IF(T .EQ. O)THEN 
Ши ту ( 
ELSEIF(T. EQ. N+K)THEN 
AT(T) = X1(T-1) 


ELSE 
АТ(Т) = КИТ) + X1(T-1) - ALPHAS(K)*XO(T-1) 
ENDIF 
60 CONTINUE 
C LOOPS TO UPDATE PREDICTION ERROR VECTORS FOR NEXT ITERATION. 
C (SHIFT X1 INTO XO AND AT INTO X1) 


00 ٣ 
ХОСТ) = X1(T) 
70 CONTINUE 
DO 80 T=0,N+K 
Х1(Т) - АТ(Т) 
80 CONTINUE 
ШЕГЕ ЕО RODD) TEEN 
KODD = KODD + 2 
ELSE 
KEVEN = KEVEN + 2 
ENDIF 
101 CONTINUE 
RETURN 
END 
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